include constants.doh

//////////////////////////////////////////////////
//// 1. Baseline March 15th Based Regressions ////
//////////////////////////////////////////////////

// See the README for instructions on how to access this data
import delimited "exposure_mvmt_regressions.csv", clear

// Make the age variables
gen age_bracket_18_34 = age_bracket == "18-24" | age_bracket == "25-34"
gen age_bracket_35_54 = age_bracket == "35-44" | age_bracket == "45-54"
gen age_bracket_55 = age_bracket == "55-64" | age_bracket == "65+"

// Scale LHS and make both positive
replace diff_share_days_no_mvmt = diff_share_days_no_mvmt*100
gen pct_reduc_avg_daily_tiles = pct_change_avg_daily_bing_tiles_*-100

// Make clean buckets
encode age_bracket, gen(age_bracket_clean)

// Make log RHS
gen log_fwc_jh = log(frnd_weighted_cases_jh_early)
gen log_fwc_100k_jh = log(frnd_weighted_cases_per_100k_jh_)

// Make the zcta/county terciles
gen top_tercile_zcta_inc = zcta_median_hh_inc_tercile == 1
gen mid_tercile_zcta_inc = zcta_median_hh_inc_tercile == 2
gen bot_tercile_zcta_inc = zcta_median_hh_inc_tercile == 3

gen top_tercile_case100k = county_cases_per_100k_jh_tercile == 1
gen mid_tercile_case100k = county_cases_per_100k_jh_tercile == 2
gen bot_tercile_case100k = county_cases_per_100k_jh_tercile == 3

// Make friend weighted measures
egen fw_median_hh_inc_grp = cut(frnd_weighted_median_hh_inc), group(100)
egen fw_pop_density_grp = cut(frnd_weighted_pop_density), group(100)
egen fw_frac_urban_grp = cut(frnd_weighted_frac_urban), group(100)

// Make interaction group for binscatter
fegen interaction_group = group(home_zcta is_female has_college age_bracket_clean has_iphone has_tablet)

// The share abroads can be NULL in rare cases.
// This forces it to 0 for these (which we believe is their true value).
// We also scale up by 100 to make it percent.
destring share_frnds_*, replace force
foreach var of varlist share_frnds_china share_frnds_korea share_frnds_italy share_frnds_spain {
    replace `var' = 0 if `var' == .
    replace `var' = `var' * 100
}

// Some of these friend rank RHS can often be NULL if you don't have friends
// that meet the criteria. If there are NULL values, will be read in as string
foreach x in friends_1_25_weighted_cases_jh_e friends_26_50_weighted_cases_jh_ friends_51_75_weighted_cases_jh_ friends_76_100_weighted_cases_jh{
    local vartype: type `x'
    if substr("`vartype'",1,3)=="str" {
        destring `x', replace force
    }
}

gen log_fwc_jh_1_25 = log(friends_1_25_weighted_cases_jh_e)
gen log_fwc_jh_26_50 = log(friends_26_50_weighted_cases_jh_)
gen log_fwc_jh_51_75 = log(friends_51_75_weighted_cases_jh_)
gen log_fwc_jh_76_100 = log(friends_76_100_weighted_cases_jh)

compress

/////// Binscatters (Figures A12 and A13) ////////

// Here we use Sergio's mata classes from REGHDFE to demean
// See: help reghdfe_mata or https://github.com/sergiocorreia/reghdfe/issues/208

preserve

    // Our raw mata commands will not be as smart about missing observations as reghdfe
    drop if log_fwc_jh == . | days_of_week != "all"

    //// Share at home LHS ////

    mata: HDFE = fixed_effects("interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp")
    mata: data = HDFE.partial_out("diff_share_days_no_mvmt log_fwc_jh")
    mata: st_store(HDFE.sample, st_addvar("double", tokens("R_diff_share_days_no_mvmt R_log_fwc_jh")), data)

    // Add back the means
    egen mean_diff_share_days_no_mvmt = mean(diff_share_days_no_mvmt)
    egen mean_log_fwc_jh = mean(log_fwc_jh)

    gen Rmean_diff_share_days_no_mvmt = R_diff_share_days_no_mvmt + mean_diff_share_days_no_mvmt
    gen Rmean_log_fwc_jh = R_log_fwc_jh + mean_log_fwc_jh

    binscatter Rmean_*, ytitle(Change in Share Days Single Tile (%)) xtitle(log(Friend Weighted Cases)) reportreg n(50) line(qfit) ytitle(, size(large)) ylabel(, labsize(large)) xtitle(, size(large)) xlabel(, labsize(large)), ///
        if days_of_week == "all"

    graph export "output/figures/binscatters/share_mvmt_logfwc.pdf"

    drop mean_diff_share_days_no_mvmt
    drop mean_log_fwc_jh
    drop R_*
    drop Rmean_*

    //// Percent Change Daily tiles LHS ////

    mata: HDFE = fixed_effects("interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp")
    mata: data = HDFE.partial_out("pct_reduc_avg_daily_tiles log_fwc_jh")
    mata: st_store(HDFE.sample, st_addvar("double", tokens("R_pct_reduc_avg_daily_tiles R_log_fwc_jh")), data)

    // Add back the means
    egen mean_reduc_avg_tiles = mean(pct_reduc_avg_daily_tiles)
    egen mean_log_fwc_jh = mean(log_fwc_jh)

    gen Rmean_pct_reduc_avg_daily_tiles = R_pct_reduc_avg_daily_tiles + mean_reduc_avg_tiles
    gen Rmean_log_fwc_jh = R_log_fwc_jh + mean_log_fwc_jh

    binscatter Rmean_*, ytitle(Reduction in Avg Daily Tiles (%)) xtitle(log(Friend Weighted Cases)) reportreg n(50) line(qfit) ytitle(, size(large)) ylabel(, labsize(large)) xtitle(, size(large)) xlabel(, labsize(large)), ///
        if days_of_week == "all"

    graph export "output/figures/binscatters/pct_reduc_daily_mvmt_logfwc.pdf"

    drop mean_reduc_avg_tiles
    drop mean_log_fwc_jh
    drop R_*
    drop Rmean_*

restore


/////// Make the determinants of friend weighted cases table (Table A5) /////////

eststo clear

eststo panel1: reghdfe log_fwc_jh age_bracket_35_54 age_bracket_55 is_female has_college has_iphone has_tablet mid_tercile_zcta_inc top_tercile_zcta_inc mid_tercile_case100k top_tercile_case100k, noabsorb vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel2: reghdfe log_fwc_jh, absorb(i.home_zcta) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel3: reghdfe log_fwc_jh, absorb(i.home_zcta fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel4: reghdfe log_fwc_jh age_bracket_35_54 age_bracket_55 is_female has_college has_iphone has_tablet, absorb(i.home_zcta fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel5: reghdfe log_fwc_jh, absorb(i.interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

xml_tab panel*, save("output/tables/covac_mvmt_constant_deters_of_fwc") stats(N r2 ymean) replace below

eststo clear


/////// Make the outcomes table (Appendix A8) ///////

eststo clear

eststo panel1: reghdfe diff_share_days_no_mvmt age_bracket_35_54 age_bracket_55 is_female has_college has_iphone has_tablet mid_tercile_zcta_inc top_tercile_zcta_inc mid_tercile_case100k top_tercile_case100k, noabsorb vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel2: reghdfe diff_share_days_no_mvmt age_bracket_35_54 age_bracket_55 is_female has_college has_iphone has_tablet, absorb(i.home_zcta) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel3: reghdfe diff_share_days_no_mvmt log_fwc_jh, absorb(i.home_zcta fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel4: reghdfe diff_share_days_no_mvmt log_fwc_jh age_bracket_35_54 age_bracket_55 is_female has_college has_iphone has_tablet, absorb(i.home_zcta fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel5: reghdfe diff_share_days_no_mvmt log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel6: reghdfe diff_share_days_no_mvmt log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "weekend"
estadd ysumm

eststo panel7: reghdfe diff_share_days_no_mvmt log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "weekday"
estadd ysumm

eststo panel8: reghdfe diff_share_days_no_mvmt log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp college) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all" & has_college == 1 & total_users_college >= 100
estadd ysumm

xml_tab panel*, save("output/tables/covac_mvmt_constant_diff_share_mvmt") stats(N r2 ymean) replace below

eststo clear


/////// Same outcome table, with the Bing tiles outcome (Appendix A23) ///////

eststo panel1: reghdfe pct_reduc_avg_daily_tiles age_bracket_35_54 age_bracket_55 is_female has_college has_iphone has_tablet mid_tercile_zcta_inc top_tercile_zcta_inc mid_tercile_case100k top_tercile_case100k, noabsorb vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel2: reghdfe pct_reduc_avg_daily_tiles age_bracket_35_54 age_bracket_55 is_female has_college has_iphone has_tablet, absorb(i.home_zcta) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel3: reghdfe pct_reduc_avg_daily_tiles log_fwc_jh, absorb(i.home_zcta fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel4: reghdfe pct_reduc_avg_daily_tiles log_fwc_jh age_bracket_35_54 age_bracket_55 is_female has_college has_iphone has_tablet, absorb(i.home_zcta fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel5: reghdfe pct_reduc_avg_daily_tiles log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel6: reghdfe pct_reduc_avg_daily_tiles log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "weekend"
estadd ysumm

eststo panel7: reghdfe pct_reduc_avg_daily_tiles log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "weekday"
estadd ysumm

eststo panel8: reghdfe pct_reduc_avg_daily_tiles log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp college) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all" & has_college == 1 & total_users_college >= 100
estadd ysumm

xml_tab panel*, save("output/tables/covac_mvmt_constant_reduc_avg_daily_bing_tiles") stats(N r2 ymean) replace below

eststo clear


//////// Make the heterogeneity tables (Tables A15 and A16) ///////

// Make all the new necessary vars
gen log_fwcXage_18_34 = log_fwc_jh*age_bracket_18_34
gen log_fwcXage_35_54 = log_fwc_jh*age_bracket_35_54
gen log_fwcXage_55 = log_fwc_jh*age_bracket_55

gen log_fwcXfemale = log_fwc_jh*(is_female==1)
gen log_fwcXmale = log_fwc_jh*(is_female==0)

gen log_fwcXcollege = log_fwc_jh*(has_college==1)
gen log_fwcXno_college = log_fwc_jh*(has_college==0)

gen log_fwcXbot_inc = log_fwc_jh*bot_tercile_zcta_inc
gen log_fwcXmid_inc = log_fwc_jh*mid_tercile_zcta_inc
gen log_fwcXtop_inc = log_fwc_jh*top_tercile_zcta_inc

gen log_fwcXbot_case100k = log_fwc_jh*bot_tercile_case100k
gen log_fwcXmid_case100k = log_fwc_jh*mid_tercile_case100k
gen log_fwcXtop_case100k = log_fwc_jh*top_tercile_case100k

compress

eststo clear

eststo panel1: reghdfe diff_share_days_no_mvmt log_fwcXage_18_34 log_fwcXage_35_54 log_fwcXage_55, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel2: reghdfe diff_share_days_no_mvmt log_fwcXfemale log_fwcXmale, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel3: reghdfe diff_share_days_no_mvmt log_fwcXcollege log_fwcXno_college, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel4: reghdfe diff_share_days_no_mvmt log_fwcXbot_inc log_fwcXmid_inc log_fwcXtop_inc, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel5: reghdfe diff_share_days_no_mvmt log_fwcXbot_case100k log_fwcXmid_case100k log_fwcXtop_case100k, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel6: reghdfe diff_share_days_no_mvmt log_fwc_jh_1_25 log_fwc_jh_26_50 log_fwc_jh_51_75 log_fwc_jh_76_100, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
        if days_of_week == "all"
estadd ysumm
test log_fwc_jh_1_25 = log_fwc_jh_76_100
estadd scalar F_top_bot = r(F)
estadd scalar p_top_bot = r(p)

xml_tab panel*, save("output/tables/covac_mvmt_constant_diff_share_mvmt_heteros") stats(N r2 ymean F_top_bot p_top_bot) replace below

eststo panel1: reghdfe pct_reduc_avg_daily_tiles log_fwcXage_18_34 log_fwcXage_35_54 log_fwcXage_55, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel2: reghdfe pct_reduc_avg_daily_tiles log_fwcXfemale log_fwcXmale, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel3: reghdfe pct_reduc_avg_daily_tiles log_fwcXcollege log_fwcXno_college, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel4: reghdfe pct_reduc_avg_daily_tiles log_fwcXbot_inc log_fwcXmid_inc log_fwcXtop_inc, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel5: reghdfe pct_reduc_avg_daily_tiles log_fwcXbot_case100k log_fwcXmid_case100k log_fwcXtop_case100k, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel6: reghdfe pct_reduc_avg_daily_tiles log_fwc_jh_1_25 log_fwc_jh_26_50 log_fwc_jh_51_75 log_fwc_jh_76_100, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
        if days_of_week == "all"
estadd ysumm
test log_fwc_jh_1_25 = log_fwc_jh_76_100
estadd scalar F_top_bot = r(F)
estadd scalar p_top_bot = r(p)

xml_tab panel*, save("output/tables/covac_mvmt_constant_reduc_avg_daily_bing_tiles_heteros") stats(N r2 ymean F_top_bot p_top_bot) replace below

eststo clear

save "data/intermediate/exposure_mvmt_regressions"

///// Appendix robustness table with international and excluding 100 miles (Table A9) /////

/// First clean the excluding 100 miles data

// See the README for instructions on how to access this data
import delimited "exposure_mvmt_regressions_100miles.csv", clear

// Scale LHS and make both positive
replace diff_share_days_no_mvmt = diff_share_days_no_mvmt*100
gen pct_reduc_avg_daily_tiles = pct_change_avg_daily_bing_tiles_*-100

// Make clean buckets
encode age_bracket, gen(age_bracket_clean)

// Make log RHS
gen log_fwc_jh = log(frnd_weighted_cases_jh_early)

// Make friend weighted measures
egen fw_median_hh_inc_grp = cut(frnd_weighted_median_hh_inc), group(100)
egen fw_pop_density_grp = cut(frnd_weighted_pop_density), group(100)
egen fw_frac_urban_grp = cut(frnd_weighted_frac_urban), group(100)

// Make interaction group
fegen interaction_group = group(home_zcta is_female has_college age_bracket_clean has_iphone has_tablet)

save "data/intermediate/exposure_mvmt_regressions_100miles"

///// Then create the table

// First the 0 distance restriction
use "data/intermediate/exposure_mvmt_regressions", clear

eststo clear

eststo panel1: reghdfe diff_share_days_no_mvmt log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

// Then the 100 distance restriction
use "data/intermediate/exposure_mvmt_regressions_100miles", clear

eststo panel2: reghdfe diff_share_days_no_mvmt log_fwc_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

// Now back to the 0 distance restriction
use "data/intermediate/exposure_mvmt_regressions", clear

eststo panel3: reghdfe diff_share_days_no_mvmt log_fwc_100k_jh, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel4: reghdfe diff_share_days_no_mvmt log_fwc_jh share_frnds_china, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel5: reghdfe diff_share_days_no_mvmt log_fwc_jh share_frnds_korea, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel6: reghdfe diff_share_days_no_mvmt log_fwc_jh share_frnds_italy, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel7: reghdfe diff_share_days_no_mvmt log_fwc_jh share_frnds_spain, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

eststo panel8: reghdfe diff_share_days_no_mvmt log_fwc_jh share_frnds_china share_frnds_korea share_frnds_italy share_frnds_spain, absorb(interaction_group fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp) vce(cluster i.home_zcta) tol(1e-5), ///
    if days_of_week == "all"
estadd ysumm

xml_tab panel*, save("output/tables/covac_mvmt_constant_appendix_robustness") stats(N r2 ymean) replace below


///////////////////////////////////////////
//// 2. Main specifications for D-in-D ////
//////////////////////////////////////////

///// Figure 2, Appendix Figures A5, A6, and A7 /////


local days_of_weeks "all weekend weekday"
local days_of_weeks_count : word count `days_of_weeks'

forvalues i = 1/`days_of_weeks_count' {

    local curr_dow : word `i' of `days_of_weeks'

    // See the README for instructions on how to access this data
    import delimited "exposure_mvmt_regressions_d_in_d_`curr_dow'.csv", clear

    // Do the encodes
    encode age_bracket, gen(age_bracket_clean)

    // Make the zcta/county terciles
    // These aren't used in the baseline, but are used in the heterogeneities
    gen top_tercile_zcta_inc = zcta_median_hh_inc_tercile == 1
    gen mid_tercile_zcta_inc = zcta_median_hh_inc_tercile == 2
    gen bot_tercile_zcta_inc = zcta_median_hh_inc_tercile == 3

    gen top_tercile_case100k = county_cases_per_100k_jh_tercile == 1
    gen mid_tercile_case100k = county_cases_per_100k_jh_tercile == 2
    gen bot_tercile_case100k = county_cases_per_100k_jh_tercile == 3

    // Make the friend weighted percentiles
    egen fw_median_hh_inc_grp = cut(frnd_weighted_median_hh_inc), group(100)
    egen fw_pop_density_grp = cut(frnd_weighted_pop_density), group(100)
    egen fw_frac_urban_grp = cut(frnd_weighted_frac_urban), group(100)

    // Make dummies for all but week6. This forces week6 to be the one omitted.
    foreach week of numlist `min_week'/`max_week' {
        gen t50_week`week' = week_of_2020 == `week' & pct_rank_fwc_zcta_early > 0.5
    }

    // Destring college so we can use it as a categorical control
    destring college, replace force

    save "data/intermediate/exposure_mvmt_regressions_d_in_d_`curr_dow'"

    // Drop a bunch of variables to make these regressions run faster
    keep share_days_no_mvmt avg_daily_bing_tiles_visited userid total_users_college college `t50_vars' week_of_2020 home_zcta fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp age_bracket_clean is_female has_college has_tablet has_iphone

    // All the RHS variables are floats -- we can compress them to simple bytes and free up a TON of memory
    compress

    //// Share days no movement ////

    reghdfe share_days_no_mvmt `t50_vars', absorb(`baseline_FEs') vce(cluster i.home_zcta) tol(1e-5)
    output_coeffs "`t50_vars'" "figures/coeffs_diff_0dist_`curr_dow'_share_days_no_mvmt"

    //// Average Bing tiles visited ////

    reghdfe avg_daily_bing_tiles_visited `t50_vars', absorb(`baseline_FEs') vce(cluster i.home_zcta) tol(1e-5)
    output_coeffs "`t50_vars'" "figures/coeffs_diff_0dist_`curr_dow'_avg_daily_bing_tiles_visited"

    //// Share days no movement - College only, controlling for specific college ////

    reghdfe share_days_no_mvmt `t50_vars', absorb(`college_FEs') vce(cluster i.home_zcta) tol(1e-5), ///
        if has_college == 1 & total_users_college >= 100
    output_coeffs "`t50_vars'" "figures/coeffs_diff_0dist_`curr_dow'_share_days_no_mvmt_college_only"

    //// Average Bing tiles visited - College only, controlling for specific college ////

    reghdfe avg_daily_bing_tiles_visited `t50_vars', absorb(`college_FEs') vce(cluster i.home_zcta) tol(1e-5), ///
            if has_college == 1 & total_users_college >= 100
    output_coeffs "`t50_vars'" "figures/coeffs_diff_0dist_`curr_dow'_avg_daily_bing_tiles_visited_college_only"
}


///////////////////////////////////
//// 3. D-in-D heterogeneities ////
///////////////////////////////////

///// Appendix Figures A8 and A9 //////

// The raw csv has been cleaned into this dta by step 2, above
use "data/intermediate/exposure_mvmt_regressions_d_in_d_all", clear

// Demean the LHS
mata: HDFE = fixed_effects("`baseline_FEs'")
mata: HDFE.tolerance = 1e-5

mata: data = HDFE.partial_out("share_days_no_mvmt avg_daily_bing_tiles_visited")
mata: st_store(HDFE.sample, st_addvar("double", tokens("R_share_days_no_mvmt R_avg_daily_bing_tiles")), data)

// Our LHS is now R_share_days_no_mvmt, the residualized version of share_days_no_mvmt
drop share_days_no_mvmt
drop avg_daily_bing_tiles_visited

//// Gender /////
preserve

    local t50_vars_gender = ""
    foreach w of numlist `min_week'/`max_week' {

        gen t50_week`w'Xfemale = t50_week`w' * (is_female==1)
        gen t50_week`w'Xmale = t50_week`w' * (is_female==0)

        local t50_vars_gender = "`t50_vars_gender't50_week`w'Xfemale t50_week`w'Xmale "
    }

    keep R_share_days_no_mvmt R_avg_daily_bing_tiles `t50_vars_gender' userid week_of_2020 fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp has_tablet has_iphone home_zcta age_bracket_clean is_female has_college

    compress

    reghdfe R_avg_daily_bing_tiles `t50_vars_gender', compact poolsize(15) absorb(`baseline_FEs') vce(cluster i.home_zcta) tol(1e-5)
    output_coeffs "`t50_vars_gender'" "figures/coeffs_diff_0dist_all_avg_daily_bing_tiles_visited_gender"

restore

//// Age ////
preserve

    local t50_vars_age = ""
    foreach w of numlist `min_week'/`max_week' {

        gen t50_week`w'X18_34 = t50_week`w' * (age_bracket == "18-24" | age_bracket == "25-34")
        gen t50_week`w'X35_54 = t50_week`w' * (age_bracket == "35-44" | age_bracket == "45-54")
        gen t50_week`w'X55 = t50_week`w' * (age_bracket == "55-64" | age_bracket == "65+")

        local t50_vars_age = "`t50_vars_age't50_week`w'X18_34 t50_week`w'X35_54 t50_week`w'X55 "

    }

    keep R_share_days_no_mvmt R_avg_daily_bing_tiles `t50_vars_age' userid week_of_2020 fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp has_tablet has_iphone home_zcta age_bracket_clean is_female has_college

    compress

    reghdfe R_avg_daily_bing_tiles `t50_vars_age', compact poolsize(15) absorb(`baseline_FEs') vce(cluster i.home_zcta) tol(1e-5)
    output_coeffs "`t50_vars_age'" "figures/coeffs_diff_0dist_all_avg_daily_bing_tiles_visited_age"

restore

//// Cases per 100k ////
preserve

    local t50_vars_cases = ""
    foreach w of numlist `min_week'/`max_week' {

        gen t50_week`w'Xbot_case = t50_week`w' * bot_tercile_case100k
        gen t50_week`w'Xmid_case = t50_week`w' * mid_tercile_case100k
        gen t50_week`w'Xtop_case = t50_week`w' * top_tercile_case100k

        local t50_vars_cases = "`t50_vars_cases't50_week`w'Xbot_case t50_week`w'Xmid_case t50_week`w'Xtop_case "

    }

    keep R_share_days_no_mvmt R_avg_daily_bing_tiles `t50_vars_cases' userid week_of_2020 fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp has_tablet has_iphone home_zcta age_bracket_clean is_female has_college

    compress

    reghdfe R_avg_daily_bing_tiles `t50_vars_cases', compact poolsize(15) absorb(`baseline_FEs') vce(cluster i.home_zcta) tol(1e-5)
    output_coeffs "`t50_vars_cases'" "figures/coeffs_diff_0dist_all_avg_daily_bing_tiles_visited_cases"

restore

//// Has College ////
preserve

    local t50_vars_college = ""
    foreach w of numlist `min_week'/`max_week' {

        gen t50_week`w'Xcollege = t50_week`w' * (has_college==1)
        gen t50_week`w'Xno_college = t50_week`w' * (has_college==0)

        local t50_vars_college = "`t50_vars_college't50_week`w'Xcollege t50_week`w'Xno_college "

    }

    keep R_share_days_no_mvmt R_avg_daily_bing_tiles `t50_vars_college' userid week_of_2020 fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp has_tablet has_iphone home_zcta age_bracket_clean is_female has_college

    compress

    reghdfe R_avg_daily_bing_tiles `t50_vars_college', compact poolsize(15) absorb(`baseline_FEs') vce(cluster i.home_zcta) tol(1e-5)
    output_coeffs "`t50_vars_college'" "figures/coeffs_diff_0dist_all_avg_daily_bing_tiles_visited_college"

restore

//// ZCTA Income ////
preserve

    local t50_vars_income = ""
    foreach w of numlist `min_week'/`max_week' {

        gen t50_week`w'Xbot_inc = t50_week`w' * bot_tercile_zcta_inc
        gen t50_week`w'Xmid_inc = t50_week`w' * mid_tercile_zcta_inc
        gen t50_week`w'Xtop_inc = t50_week`w' * top_tercile_zcta_inc

        local t50_vars_income = "`t50_vars_income't50_week`w'Xbot_inc t50_week`w'Xmid_inc t50_week`w'Xtop_inc "

    }

    keep R_share_days_no_mvmt R_avg_daily_bing_tiles `t50_vars_income' userid week_of_2020 fw_median_hh_inc_grp fw_pop_density_grp fw_frac_urban_grp has_tablet has_iphone home_zcta age_bracket_clean is_female has_college

    compress

    reghdfe R_avg_daily_bing_tiles `t50_vars_income', compact poolsize(15) absorb(`baseline_FEs') vce(cluster i.home_zcta) tol(1e-5)
    output_coeffs "`t50_vars_income'" "figures/coeffs_diff_0dist_all_avg_daily_bing_tiles_visited_income"

restore
